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The first function of the skin is to serve as a protective barrier against the 
environment. Its loss of integrity as a result of injury or illness may lead to a 
major disability and the first goal of healing is wound closure involving 
many biological processes for repair and tissue regeneration. In vivo wound 
healing has four phases, one of them being the migration of the healthy 
epithelium surrounding the wound in the direction of the injury in order to 
cover it. Here, we present a theoretical model of the re-epithelialization 
phase driven by chemotaxis for a circular wound. This model takes into 
account the diffusion of chemoattractants both in the wound and the neigh- 
bouring tissue, the uptake of these molecules by the surface receptors of 
epithelial cells, the migration of the neighbour epithelium, the tension and 
proliferation at the wound border. Using a simple Darcy's law for cell 
migration transforms our biological model into a free-boundary problem, 
which is analysed in the simplified circular geometry leading to explicit 
solutions for the closure and making stability analysis possible. It turns 
out that for realistic wound sizes of the order of centimetres and from 
experimental data, the re-epithelialization is always an unstable process and 
the perfect circle cannot be observed, a result confirmed by fully nonlinear 
simulations and in agreement with experimental observations. 



1. Introduction 

Adult skin is made up of three layers: the epidermis and the dermis separated by 
the basement membrane. When a profound injury occurs which destroys a part 
of the dermis, it has to be mended rapidly to restore the protective barrier [1,2]. 
In vivo wound healing is a complex repairing process orchestrated by intra- 
and intercellular pathways. To recover the integrity of the skin, this process is 
accomplished by four successive but overlapping stages: clotting, inflammation, 
re-epithelialization and remodelling. Immediately after an injury, a clot composed 
mainly by fibrin fibres and platelets is formed to plug the wound (stage 1). Sub- 
sequently in the next 2-10 days (stage 2) the clot is continuously infiltrated by 
inflammatory cells which clear up the debris and release chemical factors, such 
as vascular endothelial growth factor and transforming growth factor-beta. 
Once the chemical gradient is established, different cells are recruited to fabricate 
the new tissue (stage 3). On the one hand, neo-vasculature [3] and fibro- 
blasts (depositing collagen) are recruited in the dermis, reforming the clot into 
a granulation tissue, which nutritiously and physically support the repair of 
upper layers. On the other hand, keratinocytes migrate and proliferate at the 
edge of the wound to extend the newly formed epithelial carpet made of several 
layers of cells in the epidermis. This process is called re-epithelialization and 
it lasts for two to three weeks. At the end of stage 3, myofibroblasts trans- 
formed from fibroblasts contract and try to bring the wound edge together. 
They disappear by apoptosis in the dermis. The remodelling (stage 4) continues 
for months or even for years in order to restore the homoeostasis of the normal 
skin. However, the normal anatomical structure is not truly recovered and a 
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Figure 1. The closing process after a 5 mm diameter punch biopsy on the 
normal skin of wild-type mice (^7(i), from [14], copyright with permission) 
and a 6 mm diameter punch biopsy on an artificial corneal tissue (^(i), 
from [15], copyright with permission). Both present a wavy border as the 
wound closes. The photographs taken 3 days (for the wild-type mice skin) 
and 2 days (for the corneal tissue) after the punch are shown in a(\\) and 
b{\\), respectively. 



scar is formed from the granulation tissue. While this is 
true for large wounds for adult humans, those in human 
embryos [4] may be perfectly closed, the reasons being not 
fully understood. Interestingly, the scar formation is believed 
to be an evolutionary sacrifice in order to achieve rapid 
wound closure [5] and this is indicated by the intense and 
redundant inflammatory response [6] which mediates the 
systematic cell behaviours [7] during re-epithelialization 
through releasing chemical signals. In all cases, the final scar 
is an outcome of the interplay between cells and the micro- 
environment [8] through physical and chemical factors 
during the four stages. 

Here, we focus on re-epithelialization when keratinocytes, 
the epidermal predominant cells, detach from the basal mem- 
brane, migrate towards the wound margin and proliferate, 
driven by morphogens in the granulation tissue [9]. Very 
careful in vitro epithelium migration experiments have been 
performed recently by several groups [10-13]. Realized on 
a solid substrate, they involve an advancing monolayer 
motivated by the collective behaviour of cells through 
migration, generation of forces and reorganization of the cell 
cytoskeleton at the margin. The findings indicate a non-intui- 
tive manner of cells to achieve robustness during the repair 
where chemotaxis seems not to be decisive for the global 
cell migration in vitro. Nevertheless, owing to the abundant 
morphogenic signals released in the wound in vivo compared 
with in vitro experiments, chemotaxis dominates the collective 
behaviour of cells during the wound-healing process. Closer to 
in vivo wound healing (figure la(i)(ii)), an artificial cornea has 
been reconstituted [15] (figure lb{i){ii)) by the introduction of 
epithelial cells, fibroblasts and growth factors in a three-dimen- 
sional Petri dish. Interestingly, all these works [11,12,15] 
concern the circular geometry rather than the linear cut done 
in earlier experiments [16,17]. Although well controlled, it is 
unclear that these systems present the whole complexity 
of in vivo wound healing [14]. Nevertheless, both in vivo and 



in vitro experiments present disordered border evolution, 
which indicates a possible universal irregularity owing to 
chemotaxis (figure 1). 

In this work, we describe the re-epithelialization by consid- 
ering a circular hole in a tissue, adapting a recent dynamical 
model of chemotactic migration [18] driven by morphogenetic 
gradients. Being continuous, the model describes the tissue at 
a scale larger than the cell size but smaller than the size of 
the hole. Here, we do not distinguish cell species or the che- 
moattractant categories. The latter is continuously supplied 
by an incoming flux from the third dimension normal to the 
epithelial layer and uptaken by cell receptors. As shown in 
[18], the thinness of the moving layer transforms the three- 
dimensional process into a two-dimensional model, where 
the localized thickness variation at the border contributes 
to a tension T. In addition, we assume that a strong viscous 
friction exists between the moving layer and the substrate 
[19-23], which allows writing a Darcy's law for the average 
velocity field inside the tissue, while cell -cell interactions 
and mitosis are transformed into interface boundary con- 
ditions in the sharp interface limit. We first present the 
physical model, then the analytical treatment giving an explicit 
solution for the circular geometry and a study of its stability. 
Indeed in the case of a linear border, it has been shown that 
a long wavelength instability occurs at low velocities owing 
to a Goldstone mode (translational invariance) that surface ten- 
sion alone cannot succeed to prohibit [18]. Intuitively, one may 
expect that it disappears for small holes as the surface tension 
effect accentuates. On the contrary, as for experiments in vitro, 
larger holes may exhibit contour instabilities leading to dyna- 
mical behaviours that we aim to predict. To go beyond the 
stability analysis, the time-dependent free-boundary problem 
requires numerical methods able to solve equations on evol- 
ving domains. We choose to tackle our problem by level-set 
methods, first developed in the 1980s by Osher & Sethian 
[24], which are able to track moving interfaces and topological 
changes automatically and have been successfully applied to 
problems, such as liquid -gas interactions [25], image proces- 
sing [26] and tumour growth [27-30]. Owing to the nature of 
our theoretical problem, we select the methodology developed 
for tumour growth [27-30]. In the following, we first present 
the model, introducing the most relevant physical parameters, 
then the analysis for weak amplitude instability, and finally 
numerical methods and simulations for the ultimate state of 
the closure. 



2. The model 

We consider a circular hole il in a two-dimensional cellular 
population of density p immersed in a morphogenetic environ- 
ment (figure 2). Cells migrate up the morphogen gradient, 
the chemotactic flux being proportional to the concentration 
gradient of morphogens T = AcVc, with being a mobility 
constant. The morphogen repartition c(r, t) inside and outside 
the hole il is alimented through sources coming from below. 
Mitosis happens only at the border [15,31] and the closure is 
mostly achieved by the migration. The mass balance for the 
cell population is then 

p + Vpv= 7P + V - J". (2.1) 

Neglecting volumetric growth ( y = 0) in the tissue except at the 
periphery, the cell density is constant and p = po. Equation (2.1) 
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Figure 2. The schematic description: O is the wound (clot granulation 
tissue) where chemoattractants are released and il^ is the moving epithelial 
continuum. The hydrostatic pressure gradient in is driven by chemotaxis. 
The proliferation is constrained at the border 9/2 (contour in black and 
zoomed in by grey). 



simplifies to V(pv — Ac Vc) = 0 and gives that the nor- 
mal velocity of the front is directly proportional to the normal 
concentration gradient. At extremely low velocities, the cell 
migration satisfies a Darcy's law v = — Mp Vp, Mp being a por- 
osity coefficient equal to the square of the epithelium height 
divided by a friction coefficient. As shown in [18], this law is 
deduced in tissues when friction between phases or friction 
with a substrate balances the hydrostatic part of the elastic 
stress acting on the cells [18,20]. The wound perturbs the homo- 
eostatic state of the surrounding and a source of morphogens 
will try to restore an optimal equilibrium value in the aperture 
Cq. Taking this concentration as unit for the morphogen concen- 
tration (giving c = c/cq), Tc the uptake time as time units, Lg as 
length unit (with = Dq.Tc, Dg being the diffusion coefficient 
inside the epithelium), we get 



1 Ac, 



if r < 
if r > i^, 



(2.2) 



with the index h or e referring to the inner (hole) or outer 
(epithelium) domain. 8 = Dh/ Dg represents the ratio between 
the diffusion coefficient in the hole and the tissue Dg, it is 
bigger than 1, and a gives the strength of the transverse flux 
which maintains the morphogen level inside the hole. Owing 
to the relative slowness of the cell migration, we neglect the 
time dependence in equation (2.2). Taking Dg/Mp as pressure 
unit simplifies Darcy's law into v = —Vp, where for simplicity 
we keep the same notation for p and v. From mass balance 
equation (2.1), we get the following Laplace equation coupling 
the unknown pressure p to the chemoattractant concentration c: 



^{p + Ac) = 0 withA : 



^ Dg • 



(2.3) 



So, we get p = — Ac + (/), where (/> is a holomorphic function 
which satisfies A<p = 0. Finally, we must pay special attention 
to the interface boundary conditions. For equation (2.2), they 
concern the concentration continuity and the flux discontinuity 
owing to the morphogen consumption by mitosis at the border 
Oil (for demonstration, see [18]): 



and 5NVch = NVCg + / zCh- 



(2.4) 



where N is the outer normal, r'2 is the uptake rate con- 
strained at the border discussed below together with the 



mitosis rates Fi and r'^ at the border. Capillarity fixes the 
pressure jump at the interface 

1 



7^' 



(2.5) 



which gives the free-boundary condition considering the 
geometric effect where 7Z is the local radius of curvature. 
71 is equal to the radius if the geometry is a circle. However, 
it considers the local effect at small length scale when the 
boundary is perturbed. At the same time, the velocity 
dR/ dt of the interface differs from that of the epithelium v 
by the cell proliferation at the border and we get 

Vint -N.^-N.v- Tie (2.6) 



where Fi is the mitosis rate and a the capillary number related 
to the tension T at the interface so cr= (T ■ Mp)/ (dI^^tI^^). 
Wound heaHng in vivo is dominated by two-dimensional 
migration and the morphogens are probably not the nutrients, 
so F2 vanishes and /\c must be replaced by F'^ if we also con- 
sider mitosis without the quantitative effect of the nutrient 
concentrations. The capillary number a may involve the activi- 
ties of actin cable (bundled microfilaments) in the cells. During 
embryonic wound healing, the leading edge cells can coordi- 
nate their actin cables globally which generate an effect on 
the macroscopic level. However, this coordinated behaviour 
is lost in adults, so the effect of actin cables may act only locally. 
This resembles a in our model which is not effective for large 
wounds. Considering the set of equations (2.2) and (2.3) in 
conjunction with the set of boundary conditions (2.4), (2.5) 
and (2.6) shows that chemotatic-driven migration is indeed a 
free-boundary problem [32,33] involving several parameters 
to represent the biological complexity and the study of 
simple cases, for example the circular closure, may help in 
understanding their corresponding roles. Thus, we present 
first the analytical results for a hole remaining circular at all 
times, then its stability. 



3. The results 

3.1. Regular circular closure 

In the quasi-static approximation, equation (2.2) can be 
solved analytically and gives 

if r > i^, 
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Jo (resp. Kq) being the modified Bessel function of zero order, 
regular at r = 0 (resp. r ^ 00). Equation (3.1) takes into 
account the continuity at the interface, A being given by the 
flux continuity and reads 

{c^8)JCi{R) 



A- 



{a8)K:i{R)+Xi{aR)]' 



(3.2) 



with the following definition for Xy(X) = Jy_i(X)/Jy(X) 
(ratio of two successive Bessel functions) and the equiva- 
lent for /Cy(X) = K^_i{X)/Kj{X). The pressure Pq inside the 
epithelium becomes 

'li-)-Wi)-r 
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(3.3) 
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Figure 3. The closing velocity /? as a function of/?, varying M 5, (^) a and (c) (? = kJA. For large /?, the velocity is constant. The velocity converges to 0 when /? 
goes to 0 except that the effect of yijn(>0) is considered (c). 
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We use the Laplace law to fix the unknown degree of free- 
dom. The holomorphic function (/>, proportional to log(r), 
represents a possible driving force which appears at the inter- 
face (i.e. the so-called kenotaxis defined in [13]). Indeed, even in 
the absence of morphogens, in vitro epithelium migration on 
solid substrate is observed, perhaps owing to the reorganiz- 
ation of the cellular cytoskeleton at the wound margin. Note 
that the model remains valid if we discard the chemotaxis pro- 
vided the definition of length and time units is modified. 
Furthermore, Cq and Pq are dependent on time via RJ^). For 
simplicity, we drop the time dependence of R. The velocity R 
of the closure deduced from equation (3.5) is then 

Kx(R) R ■ 



R 



(3.4) 



The closing velocity is constant for large holes R = —a8{A-\- 
ri)/(l + a8) and is exponentially small R - e-io^'s{A+rr)/2)t 
when the radius R becomes tiny, if one restricts on chemo- 
tatic migration and proliferation. So the radius of a large 
hole begins with a linear decrease in time but the total 
closure will take an infinite time for a complete achieve- 
ment. For small holes, chemotaxis becomes subdominant 
compared with kenotaxis Ai^ at the interface which controls 
the closure dynamics and the radius satisfies a diffusive 
law in t^^^ [12]. Figure 3 shows A as a function of R for 
different values of a, 8 and Ai^/ A and §4 contains a 
discussion of parameters. 

3.2. Loss of circularity 

However, equation (3.4) is valid only if the circular contour is 
maintained. It is why we perform a linear stability analysis 
for large wounds {Ai^^ 0) assuming a small harmonic pertur- 
bation for the radius as Rb{t) = R{t)[l ^ e ■ e^-^cos{n0)] 
inducing variations on the pressure P and concentration 
field C of the same order e as follows: 



C(r, 0, t) 
and P(r, 0, t) 



:Po(r, f) + f ^e^«W(n^) 



(3.5) 



where the subscript / indicates either a quantity relative to 
the hole (h, 0 < r < R) or to the epidermis (e, R < r < oo). 
Although all perturbative quantities depend on the selected 
mode n, the linear perturbation analysis treats these modes 
independently. So we drop the n index and calculate the 



perturbative concentration fields C/(r) from equation (2.3) 
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with = v/ 1 + il„ , 



while the perturbed pressure field is given by 
(r) = BR^ /r" — Acq (r) where we take into account the har- 
monic modes of the holomorphic function (/> with B fixed 
by the Laplace law B = RR + Ax^ — cr{rp- — 1)/R. Owing to 
the weakness of e, our system of equations, once linearized, 
can be solved analytically and the growth rate of the mode 
n reads 



n(n 



1) n 



■ c^RIC„{a2R)], 



(3.7) 



where 



aXijaR) + JCijR) + {8 - l){aiXn{aiR) - n/R) 
{8 - l)n/R - {8aiIn{aiR) + a2/Cn(c^2i^)) 



Equation (3.7) represents an implicit relationship for solved 
by iterative techniques and positive values indicate the modes 
responsible for the destabilization of the circular border as the 
migration proceeds. The results are presented for different 
values of A, a and 8, where il„ is displayed in figure 4. See §4 
for the discussion of parameters. The results indicate an instabil- 
ity leading to a deviation from a circle in short time, for n up to a 
critical mode (fixed by the capillarity). It is why numerical 
simulations are necessary to go beyond the linear analysis and 
fully consider the nonlinearities. 



3.3. Full dynamics and numerical methods 

We discretize c in equation (2.2) and pressure p in equation 
(2.3) on a Cartesian mesh in space and implicitly in time, 
using a nonlinear adaptive Gauss -Seidel iterative method 
[29,30]. At the domain boundary, we impose both vanishing 
values of c and normal gradient of p. Then the noise is added 
on C = c -\- X with |;^| < 1. Beginning with a perfect circle 
with R = 90, the wound closing is tracked by the level-set 
method developed in [27,28] where a scalar function (P 
with |V0| = 1 describes the wound (^>< 0), the epithelium 
(0 > 0) and the interface (0=0). The normal and curvature 
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Figure 4. Growth rate il„ as a function of wavenumber n varying [a) the mobility coefficient A, [b) the strength of morphogen a and (c) the diffusion coefficient 
ratio 8 = Dh/Z^e between the wound and the epithelium. 




Figure 5. The simulation of wound closing. Both begin with /? = 90 (9 mm) 
and 8 = 2. {a)(\) The concentration field and (ii) the overall time evolution 
with A = 1 and a = 1; (^)(i) the concentration field and (ii) the overall time 
evolution with A = 10 and a = 0.1. Further explanation is given in the text. 



in equation (2.5) are calculated by standard differential geo- 
metry: V0/|V0| and V • (VcP/|VcP|). 0 is updated by the 
relation 0t + V^ext|VcP| = 0, where Vgxt is the constant exten- 
sion of Vint on the interface from equation (2.5) following 
the gradient VcP. Note that this method simulates a perfect 
circular closing without noise. According to our results, the 
wound deviates from a perfect circle soon after the growth 
as shown in figures 5 and 6 and in agreement with the stab- 
ility analysis. 



4, Discussion and conclusion 

There are several independent parameters in this model. We 
can fix some of them with the published experimental data 
(table 1). Our velocity unit is compatible with the closing vel- 
ocity of the cornea (3 days for a hole of radius approximately 
3 mm [15]) giving A^l for this experiment according to 
equation (3.2). In figure 4, we vary the parameters a and 5, 
which are more difficult to estimate, and also A. Owing to 
the large size of wounds in practice, the linear stability analy- 
sis gives always an instability and this conclusion is robust 




Figure 6. The simulation of the late-stage wound closing with large defor- 
mation with 8=2 (top, left) and 5 = 1 (bottom, left). Right: The MDCK 
cell monolayer with a wound made by a polydimethylsiloxane cylindrical 
pillar after 10 (up) and 15 (bottom) hours [12]. The initial diameter of 
the wound in the experiment is 0.5 mm. The nonlinear border (red 
arrows) at late stage merges and leads to the pinch-off of wound area 
(orange arrows) in both the simulation and the experiment. 

to parameter changes. This finding is confirmed by fully non- 
linear simulations under the same parameter range (figures 5 
and 6). In the simulations, large A (approx. 10, total time 
approx. 60) contributes to faster closing compared with 
small A (approx. 1, total time approx. 600). A typical exam- 
ple of the closing process with a = l and 5 = 2 is shown in 
figure 5a. The advancing interface becomes wavy as the 
wound heals; however, when the wound becomes small, 
surface tension re-stabilizes the wound to a potato shape 
consistent with the linear stability analysis. This surface ten- 
sion may also pinch wounds off to smaller pieces as shown 
in figure 5h{ii), figure 6 and in [12]. This event requires at first 
a more irregular closing (red arrows in figure 6), when the con- 
centration of chemoattractant is less homogeneous given 
inadequate transverse flux {a = 0.1). The pinch-off can 
happen both at the intermediate or the end stage of the closing 
(yellow arrows), consistent with observations at the late stage 
of wound healing on a monolayer of MDCK cells (figure 6, 
right), suggesting the same behaviour in more complex 
wound-healing processes. 

In this work, we have shown theoretically that wound 
re-epithelialization gives a border instability under clinically 
realistic parameters. Driven by chemotaxis, our model does 
not introduce other cell activities which close the wound at 
the ultimate stage. However, we conjecture that this chemo- 
tactic instability at the border of the wound may affect the 
quality of the final repair. During embryonic wound healing 
where perfect reconstruction is observed, the wound is closed 



Table 1. The parameter table. 



physical parameter 


value 


diffusion coefficient 


1 (Jim s Lj4J 


uptake time Tc 


-)AAA c ['3A^ 

zUUU S Lj4J 


surface tension T 


lU IM m [IzJ 


friction coefficient 


lU IM (xm LIzJ 


the length unit 


50 |jLm, calculated 


the velocity unit 


0.025 iJim s~\ calculated 


the capillary number cr 


0.1, calculated 


the chemoattractant strength a 


0.1-10, estimated 


the diffusion coefficient ratio 8 


0.1-10, estimated 


the cell velocity A 


1-10, estimated 



dimension which we consider as 'normar. Also, as the 
wound tissue undergoes remodelling during the healing, 
the mechanical parameters cannot be considered as constants 
and the relaxation of the tissue probably evolves on a larger 
time scale than the re-epithelialization. Indeed, these relax- 
ation processes should be also considered in predicting the 
quality of the final repair. Before that, the irregularity of the 
wound border driven by chemotaxis should be considered 
as an ingredient. 
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by a 'purse string' which represents a coordination of leading 
edge cells facilitated by the actin cable [35]. This mechanism 
is lost in the adult skin, where the wound is closed by the 
crawling of cells (several layers) on the granulation tissue 
[36]. Compared with a static substrate in vitro, the granula- 
tion tissue undergoes contraction at the end of the 
re-epithelialization by myofibroblasts [37]. Indeed, this pur- 
pose is to bring the edge together which resembles the 
'purse string', but the contraction by those cells needs to be 
compatible with the synchronous material constitution, 
which is mechanically and systematically a challenge [7]. 

At the end, we discuss our model in the context of the 
skin wound healing in vivo. As the wound goes deeply into 
the dermis which is the case for most wounds, the two 
layers behave differently during re-epithelialization. The 
depth of a wound contributes to the thickness in the lower 
layer where the chemoattractant is released. Considering 
the re-epithelialization where only several layers of cells 
become motile above the lower stacks where the chemo- 
attractant is released, the depth of the wound contributes 
mathematically to the transverse chemoattractant flux 
strength a. Besides a, which may involve the depth of the 
wound, the geometry, the size of the wound as well as 
the diffusion coefficient ratio 8 are considered. Comparing 
with experiments in vivo, the human wound-healing process 
may be different from that in experimental mice owing to 
different skin biomechanical properties. However, if one 
intends to consider the problem in vivo accounting for the 
bioelasticity, the challenge does not only come from different 
components of the two layers, but also the fact that the junc- 
tion between the dermal and epidermal layer is highly 
disordered even for the intact skin [38], as shown in our pre- 
vious work. A regular pattern should be defined in the third 



Appendix A, Numerical implementation 

The simulation is done on a 200 x 200 Cartesian lattice with a 
grid size of 10 |jLm. This free-boundary problem is easy to be 
implemented with the level-set method owing to its depen- 
dence on the curvature in solving the pressure in equation 
(2.5). The full method is adapted from [28,29] and has an 
order of accuracy of more than 1.5. While the boundary can 
be implicitly given by the zero contour of the level-set 
function (p, the curvature can also be readily calculated. At 
each time step, the domain from (It to (2t+At is updated 
as follows: 

— Solve equation (2.2) with the fixed domain fit in the 
steady state with the free-boundary condition equation 
(2.4). At the boundary of the computational lattice, the 
Neumann boundary condition is imposed. 

— Update the solution of concentration C by C = c -\- x 
pointwisely, where x is randomly generated out of uni- 
form distribution between —\Xc\ arid l^'d • lAfd = ^-^^ is 
used for the presented simulations. 

— Solve equation (2.3) with known C under the interfacial 
boundary condition from the first part of equation (2.5), 
where the local curvature l/TZ is calculated from the 
level-set function by V • (V(P/|V(P|). The Neumann 
boundary condition is imposed at the boundary of the 
computational lattice. 

— Calculate the velocity V of the moving boundary from the 
second part of equation (2.4), where N is calculated from 
the level-set function by V0/|V0|. 

— Find the appropriate given by the CFL condition by 
At < Ax/4max|y| and update the domain Of to flt^M- 

— Start with the newly updated domain ilf+Af and repeat 
the last five steps. 
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